Grazing Tradeoff Analysis

Install packages

require(Rmisc) # Package for doing data summaries. We'll use its summarySE function later to get means and SEs of various datasets for plotting.

## Loading required package: Rmisc

## Loading required package: lattice

## Loading required package: plyr

require(sciplot) # The bargraph.CI function in this package makes barplots w/ error bars a breeze!

## Loading required package: sciplot

library(Rmisc)
library(sciplot)
library(plyr)

library(stringr) # for extracting numeric values

require(segmented) # for segmented line fitsß

## Loading required package: segmented

## Loading required package: MASS

## Loading required package: nlme

Import data

graz.data <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/AGCamara/Grazing_all_Och_MH_MD - ExptData_v2.csv")
#graz.data<- graz.data[graz.data$OldData!= "Y",]

# From Barbaglia et al:
Barbaglia.data <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/AGCamara/BarbagliaData/OchTradeoffData.csv")

Plotting details

# Rainbow of colours for the eight strains (as in Barbaglia et al 2024)
# Colours from: https://github.com/BlakeRMills/MetBrewer - Hokusai 1 + a bit of Hokusai 2
col.a <- rgb(96/255,26/255,16/255)
col.b <- rgb(183/255,51/255,48/255)
col.c <- rgb(234/255,95/255,75/255)
col.d <- rgb(233/255,120/255,47/255)
col.e <- rgb(243/255,186/255,81/255)
col.f <- rgb(120/255,171/255,125/255)
col.g <- rgb(0/255,93/255,148/255)
col.h <- rgb(0/255,56/255,98/255)
col.j <- rgb(64/255, 0/255, 128/255)


# Define the colors
col.a <- rgb(96/255,26/255,16/255)
col.b <- rgb(183/255,51/255,48/255)
col.c <- rgb(234/255,95/255,75/255)
col.d <- rgb(233/255,120/255,47/255)
col.e <- rgb(243/255,186/255,81/255)
col.f <- rgb(120/255,171/255,125/255)
col.g <- rgb(0/255,93/255,148/255)
col.h <- rgb(0/255,56/255,98/255)
col.j <- rgb(64/255, 0/255, 128/255)


# Rainbow of colours with reduced opacity for backgrounds of plots
alpha.bg <- .5
col.a.bg <- rgb(96/255,26/255,16/255,alpha.bg)
col.b.bg <- rgb(183/255,51/255,48/255,alpha.bg)
col.c.bg <- rgb(234/255,95/255,75/255,alpha.bg)
col.d.bg <- rgb(233/255,120/255,47/255,alpha.bg)
col.e.bg <- rgb(243/255,186/255,81/255,alpha.bg)
col.f.bg <- rgb(120/255,171/255,125/255,alpha.bg)
col.g.bg <- rgb(0/255,93/255,148/255,alpha.bg)
col.h.bg <- rgb(0/255,56/255,98/255,alpha.bg)
col.j.bg <- rgb(64/255, 0/255, 128/255, alpha.bg)

# Assembling into a vector of colors
col.strains <- c(col.a,col.b,col.c,col.d,col.e,col.f,col.g,col.h, col.j)
col.strains.bg <- c(col.a.bg,col.b.bg,col.c.bg,col.d.bg,col.e.bg,col.f.bg,col.g.bg,col.h.bg,col.j.bg)

# "Mega" vector alternates saturated and unsaturated, for making barplots that show high and low bacterial treatments side by side
mega.col.strains <- rep(NaN,length(col.strains)*2)
mega.col.strains.white <- rep(NaN,length(col.strains)*2)

for(i in 1:length(col.strains)){
  mega.col.strains[2*(i-1)+1] <- col.strains.bg[i]
  mega.col.strains.white[2*(i-1)+1] <- 'white'
  mega.col.strains[2*i] <- col.strains[i]
  mega.col.strains.white[2*i] <- col.strains[i]
}


# Line types for temperatures
ltys <- c(1,2)

# Shapes for replicates
rep.pchs <- c(1,4,5)

xenic.strains <- c(584,590,1148,1150,1391,1392,1393,2951)
rice.strains <- c('584R','590R','1148R','1150R','1391R',"1392R",'1393R','2951R')

Plot population dynamics

# Make a list of all the different prey concentrations
preyconc <- levels(as.factor(graz.data$TargetOchConc))
preyconc

##  [1] "0"      "10000"  "15000"  "20000"  "25000"  "40000"  "50000"  "60000" 
##  [9] "70000"  "75000"  "90000"  "95000"  "100000" "106000" "125000"

# Dates for plotting
days <- c(1:4)
days

## [1] 1 2 3 4

# Select the columns to use for plotting
preycolumns <- c(12:15)+1
head(graz.data[,preycolumns]) # see how this subsets the data to just the prey columns?

##   D1OchPop D2OchPop D3OchPop D4OchPop
## 1      741      228      969     4883
## 2    20511    15726    24557    33501
## 3    49282    31967    57264    63883
## 4    76973    48726    77484    87685
## 5   107225    76510    89041    85538
## 6   125564   105118   100437   102322

predcolumns <- c(17:20)+1
head(graz.data[,predcolumns])

##   D1PredPop D2PredPop D3PredPop D4PredPop
## 1         0         0         0         0
## 2         0         0         0         0
## 3         0         0         0         0
## 4         0         0         0         0
## 5         0         0         0         0
## 6         0         0         0         0

# Set up colors for each prey concentration; using VanGogh3 palette from MetBrewer https://github.com/BlakeRMills/MetBrewer
ltgreen <- rgb(156/255,193/255,132/255)
dkgreen <- rgb(31/255,90/255,37/255)
ctrlcol <- 'gray80'
colfunc <- colorRampPalette(c(ltgreen,dkgreen))
col2use <- colfunc(length(preyconc))
col2use

##  [1] "#9CC184" "#93B97D" "#8AB276" "#81AA6F" "#78A368" "#6F9C62" "#66945B"
##  [8] "#5D8D54" "#54864D" "#4B7E46" "#427740" "#397039" "#306832" "#27612B"
## [15] "#1F5A25"

# Set up line types for different predation treatments
predtreats <- levels(as.factor(graz.data$Pred))
predtreats 

## [1] "Control"  "O.marina"

ltys <- c(1:length(predtreats))

# Set up shapes for different replicates
pchs <- c(21:25)
reps <- unique(graz.data$Rep)

# Create a list of unique experiments
strains <- sort(unique(graz.data$Strain))

Calculating growth rates

# Add new columns to our data frame that we will use for prey and predator growth rates and y-intercepts. We'll consider both 48 and 72hr windows for the prey.
graz.data$prey.mu.48 <- NaN
graz.data$prey.yint.48 <- NaN
graz.data$prey.R2.48 <- NaN

graz.data$prey.mu.72 <- NaN
graz.data$prey.yint.72 <- NaN
graz.data$prey.R2.72 <- NaN

graz.data$pred.mu.72 <- NaN
graz.data$pred.yint.72 <- NaN
graz.data$pred.R2.72 <- NaN

# Fit the growth curves line-by-line, using lm()

for( i in 1:dim(graz.data)[1] ){ # for each row in the data table
  if(graz.data[i,]$AssayTemp %in% c(24)){ # excluding 30*C data
  
  if(graz.data[i,]$TargetOchConc!=0){ # Only fit prey if we expect positive densities
  # Fit all prey
  lm1 <- lm(log(t(graz.data[i,preycolumns])+1)~days)
  graz.data$prey.mu.72[i] <- summary(lm1)$coefficients[2,1]
  graz.data$prey.yint.72[i] <- summary(lm1)$coefficients[1,1]
  graz.data$prey.R2.72[i] <- summary(lm1)$r.squared
  
  lm1 <- lm(log(t(graz.data[i,preycolumns[1:3]])+1)~days[1:3])
  graz.data$prey.mu.48[i] <- summary(lm1)$coefficients[2,1]
  graz.data$prey.yint.48[i] <- summary(lm1)$coefficients[1,1]
  graz.data$prey.R2.48[i] <- summary(lm1)$r.squared
  
  }
  
  if( graz.data[i,]$Pred!='Control' ){ # Only fit predators if we expect non-zero densities
  # Fit predator
    if(!is.na(graz.data[i,]$D1PredPop)){
  lm1 <- lm(log(t(graz.data[i,predcolumns]))~days)
  graz.data$pred.mu.72[i] <- summary(lm1)$coefficients[2,1]
  graz.data$pred.yint.72[i] <- summary(lm1)$coefficients[1,1]
  graz.data$pred.R2.72[i] <- summary(lm1)$r.squared
  }}
  }
}

Let’s look at the differences in our 48 vs 72hr curve fits and check our work by adding our growth curves to our plots of population dynamics. This code creates one set of plots for every Ochromonas strain type; for strains for which multiple experiments were run, the data are broken down by experiment.

#par(mar=c(4,5,1,1),mfrow=c(4,2))


for(strainctr in 1:length(strains)){

  #if(rem(strainctr,2)==1){
    par(mar=c(4,5,1,1),mfrow=c(3,2))

 # }

StrainChoice <- strains[strainctr]
AssayLightChoice <- 100
AssayTempChoice <- 24

# Subset the dataframe to the focal Ochromonas strain
dat5 <- graz.data[graz.data$Strain==StrainChoice & graz.data$AssayLight==AssayLightChoice&graz.data$AssayTemp==AssayTempChoice,]
if(strainctr==3){ dat5 <- dat5[dat5$OldData!='Y',]}

for(datectr in 1:length(unique(dat5$D1Date))){
  dat1 <- dat5[dat5$D1Date==unique(dat5$D1Date)[datectr],]



# Create a vector of prey concentrations for this Ochromonas strain
preyconc <- levels(as.factor(dat1$TargetOchConc))

colfunc <- colorRampPalette(c('gray90',col.strains[ceiling(strainctr/2)]))
col2use <- colfunc(length(preyconc))


dat2 <- dat1[dat1$TargetOchConc!=0,]


plot(days,dat2[1,preycolumns]+1,las=1,xlab='Time (days)',ylab='',log='y',ylim=c((min(dat2[,preycolumns]+1,na.rm=T)),(max(dat2[,preycolumns]+1,na.rm=T))),type='n',main=paste('CCMP',StrainChoice, ', ', dat2$D1Date[1]))
mtext('Population size (cells/mL)',side=2,line=4)

for(i in 2:length(preyconc)){
  subdat <- dat2[dat2$TargetOchConc==preyconc[i],]
  for(j in 1:dim(subdat)[1]){
    points(days,subdat[j,preycolumns],col=col2use[i],pch=21,bg=col2use[i])
    if(subdat[j,]$Pred=='Control'){
      # 48hr growth curve
      xcoords <- seq(from = min(days), to = max(days[3]), length.out = 100)
      ycoords <- exp(subdat[j,]$prey.yint.48)*exp(subdat[j,]$prey.mu.48*xcoords)
      lines(xcoords,ycoords,col=col2use[i])
      lines(xcoords,ycoords,col='red')
      
      # 72hr growth curve
      xcoords <- seq(from = min(days), to = max(days), length.out = 100)
      ycoords <- exp(subdat[j,]$prey.yint.72)*exp(subdat[j,]$prey.mu.72*xcoords)
      lines(xcoords,ycoords,col=col2use[i])
      
    }      else{
      # 48hr growth curve
      xcoords <- seq(from = min(days), to = max(days[3]), length.out = 100)
      ycoords <- exp(subdat[j,]$prey.yint.48)*exp(subdat[j,]$prey.mu.48*xcoords)
      lines(xcoords,ycoords,col=col2use[i],lty=2)
      lines(xcoords,ycoords,col='red',lty=2)
      
      # 72hr growth curve
      xcoords <- seq(from = min(days), to = max(days), length.out = 100)
      ycoords <- exp(subdat[j,]$prey.yint.72)*exp(subdat[j,]$prey.mu.72*xcoords)
      lines(xcoords,ycoords,col=col2use[i],lty=2)
    }
  }
}


dat2 <- dat1[dat1$Pred!='Control',]

plot(days,dat2[1,predcolumns],las=1,xlab='Time (days)',ylab='Population size (cells/mL)',log='y',ylim=c((min(dat2[,predcolumns],na.rm=T)),(max(dat2[,predcolumns],na.rm=T))),type='n',main='Predator')
#mtext('Population size (cells/mL)',side=2,line=4)

for(i in 1:length(preyconc)){
  subdat <- dat2[dat2$TargetOchConc==preyconc[i],]
  for(j in 1:dim(subdat)[1]){
    points(days,subdat[j,predcolumns],col=col2use[i],pch=21,bg=col2use[i])
      xcoords <- seq(from = min(days), to = max(days), length.out = 100)
      ycoords <- exp(subdat[j,]$pred.yint.72)*exp(subdat[j,]$pred.mu.72*xcoords)
      lines(xcoords,ycoords,col=col2use[i])
      
  }
}
}
}

Compute grazing rates

Our goal is to use the growth rates of the different organisms to determine the grazing rate. A core part of this is determining differences in prey growth rates between experimental and control wells. We should see that, in the presence of O. marina, our prey grow slower.

The first thing that we need to do is match up our control (predator-free) growth rates with our predator ones.

graz.data$Strain.D1Date.TargetOchConc.Rep <- paste(graz.data$Strain,graz.data$D1Date,graz.data$TargetOchConc,graz.data$Rep,sep='.')

graz.data$ctrl.prey.mu.72 <- graz.data[graz.data$Pred=='Control',]$prey.mu.72[match(graz.data$Strain.D1Date.TargetOchConc.Rep,graz.data[graz.data$Pred=='Control',]$Strain.D1Date.TargetOchConc.Rep)] # the "match" function is a really powerful tool for linking data. Here, I'm using the "target prey" variable to make sure I match the correct average growth rates to each row in my data frame.

head(graz.data) # we can inspect the data frame to see if this worked. notice how the growth rates line up?

##   OldData OldDat2 Strain Rep AssayTemp AssayLight TargetOchConc    Pred  D1Date
## 1                    584   A        24        100             0 Control 1/30/24
## 2                    584   A        24        100         25000 Control 1/30/24
## 3                    584   A        24        100         50000 Control 1/30/24
## 4                    584   A        24        100         75000 Control 1/30/24
## 5                    584   A        24        100        100000 Control 1/30/24
## 6                    584   A        24        100        125000 Control 1/30/24
##   VolOch VolK  X D1OchPop D2OchPop D3OchPop D4OchPop X.1 D1PredPop D2PredPop
## 1     NA   NA NA      741      228      969     4883  NA         0         0
## 2     NA   NA NA    20511    15726    24557    33501  NA         0         0
## 3     NA   NA NA    49282    31967    57264    63883  NA         0         0
## 4     NA   NA NA    76973    48726    77484    87685  NA         0         0
## 5     NA   NA NA   107225    76510    89041    85538  NA         0         0
## 6     NA   NA NA   125564   105118   100437   102322  NA         0         0
##   D3PredPop D4PredPop X.2 Scope Run   prey.mu.48 prey.yint.48   prey.R2.48
## 1         0         0  NA    MD   1          NaN          NaN          NaN
## 2         0         0  NA    MD   1  0.090013793     9.720203 0.1611998785
## 3         0         0  NA    MD   1  0.075055211    10.560980 0.0614850373
## 4         0         0  NA    MD   1  0.003308333    11.094400 0.0001547864
## 5         0         0  NA    MD   1 -0.092915294    11.594080 0.3021313158
## 6         0         0  NA    MD   1 -0.111641466    11.830191 0.8954060271
##    prey.mu.72 prey.yint.72 prey.R2.72 pred.mu.72 pred.yint.72 pred.R2.72
## 1         NaN          NaN        NaN        NaN          NaN        NaN
## 2  0.19174439     9.550652  0.6074022        NaN          NaN        NaN
## 3  0.13614234    10.459168  0.3344573        NaN          NaN        NaN
## 4  0.08547347    10.957459  0.1822492        NaN          NaN        NaN
## 5 -0.05262254    11.526925  0.2341172        NaN          NaN        NaN
## 6 -0.06596196    11.754058  0.6879644        NaN          NaN        NaN
##   Strain.D1Date.TargetOchConc.Rep ctrl.prey.mu.72
## 1                 584.1/30/24.0.A             NaN
## 2             584.1/30/24.25000.A      0.19174439
## 3             584.1/30/24.50000.A      0.13614234
## 4             584.1/30/24.75000.A      0.08547347
## 5            584.1/30/24.100000.A     -0.05262254
## 6            584.1/30/24.125000.A     -0.06596196

However, because the control data are rather noisy for these experiments, we’re also going to use an alternate approach in which we use the mean growth rates from the control prey treatments and the mean growth rates for the predators. This helps to eliminate the rep-by-rep noise and give more accurate estimates of grazing rates.

graz.data$Strain.D1Date.TargetOchConc <- paste(graz.data$Strain,graz.data$D1Date,graz.data$TargetOchConc,sep='.')
ctrl.mu.summ <- summarySE(data=graz.data[graz.data$Pred=='Control',],measurevar='prey.mu.72',groupvars = 'Strain.D1Date.TargetOchConc')


graz.data$ctrl.prey.mu.72.mean <- ctrl.mu.summ$prey.mu.72[match(graz.data$Strain.D1Date.TargetOchConc,ctrl.mu.summ$Strain.D1Date.TargetOchConc)] # the "match" function is a really powerful tool for linking data. Here, I'm using the "target prey" variable to make sure I match the correct average growth rates to each row in my data frame.

head(graz.data) # we can inspect the data frame to see if this worked. notice how the growth rates line up?

##   OldData OldDat2 Strain Rep AssayTemp AssayLight TargetOchConc    Pred  D1Date
## 1                    584   A        24        100             0 Control 1/30/24
## 2                    584   A        24        100         25000 Control 1/30/24
## 3                    584   A        24        100         50000 Control 1/30/24
## 4                    584   A        24        100         75000 Control 1/30/24
## 5                    584   A        24        100        100000 Control 1/30/24
## 6                    584   A        24        100        125000 Control 1/30/24
##   VolOch VolK  X D1OchPop D2OchPop D3OchPop D4OchPop X.1 D1PredPop D2PredPop
## 1     NA   NA NA      741      228      969     4883  NA         0         0
## 2     NA   NA NA    20511    15726    24557    33501  NA         0         0
## 3     NA   NA NA    49282    31967    57264    63883  NA         0         0
## 4     NA   NA NA    76973    48726    77484    87685  NA         0         0
## 5     NA   NA NA   107225    76510    89041    85538  NA         0         0
## 6     NA   NA NA   125564   105118   100437   102322  NA         0         0
##   D3PredPop D4PredPop X.2 Scope Run   prey.mu.48 prey.yint.48   prey.R2.48
## 1         0         0  NA    MD   1          NaN          NaN          NaN
## 2         0         0  NA    MD   1  0.090013793     9.720203 0.1611998785
## 3         0         0  NA    MD   1  0.075055211    10.560980 0.0614850373
## 4         0         0  NA    MD   1  0.003308333    11.094400 0.0001547864
## 5         0         0  NA    MD   1 -0.092915294    11.594080 0.3021313158
## 6         0         0  NA    MD   1 -0.111641466    11.830191 0.8954060271
##    prey.mu.72 prey.yint.72 prey.R2.72 pred.mu.72 pred.yint.72 pred.R2.72
## 1         NaN          NaN        NaN        NaN          NaN        NaN
## 2  0.19174439     9.550652  0.6074022        NaN          NaN        NaN
## 3  0.13614234    10.459168  0.3344573        NaN          NaN        NaN
## 4  0.08547347    10.957459  0.1822492        NaN          NaN        NaN
## 5 -0.05262254    11.526925  0.2341172        NaN          NaN        NaN
## 6 -0.06596196    11.754058  0.6879644        NaN          NaN        NaN
##   Strain.D1Date.TargetOchConc.Rep ctrl.prey.mu.72 Strain.D1Date.TargetOchConc
## 1                 584.1/30/24.0.A             NaN               584.1/30/24.0
## 2             584.1/30/24.25000.A      0.19174439           584.1/30/24.25000
## 3             584.1/30/24.50000.A      0.13614234           584.1/30/24.50000
## 4             584.1/30/24.75000.A      0.08547347           584.1/30/24.75000
## 5            584.1/30/24.100000.A     -0.05262254          584.1/30/24.100000
## 6            584.1/30/24.125000.A     -0.06596196          584.1/30/24.125000
##   ctrl.prey.mu.72.mean
## 1                  NaN
## 2           0.17969735
## 3           0.15188058
## 4           0.05697058
## 5           0.02953099
## 6          -0.01727983

From here, we can compute our grazing rates as the difference between the control and +predator growth rates.

# Rep-by-rep
graz.data$prey.g.72 <- graz.data$ctrl.prey.mu.72-graz.data$prey.mu.72
graz.data$prey.g.72 <- graz.data$ctrl.prey.mu.72-graz.data$prey.mu.72

# Mean
graz.data$prey.g.72.mean <- graz.data$ctrl.prey.mu.72.mean-graz.data$prey.mu.72
graz.data$prey.g.72.mean <- graz.data$ctrl.prey.mu.72.mean-graz.data$prey.mu.72

Computing the per-predator ingestion rate is a little trickier. Ask someone to show you the calculus on the whiteboard, but essentially what we are doing is using the growth rates of both predators and prey to compute mean populations of both, and then computing the per-predator ingestion rate from there. This method follows Jeong & Latz (1994) (https://www.int-res.com/articles/meps/106/m106p173.pdf)

graz.data$prey.popn.24 <- graz.data$D1OchPop*(exp(graz.data$prey.mu.48)-1)/(graz.data$prey.mu.48) # mean prey population
graz.data[graz.data$TargetOchConc==0,]$prey.popn.24 <- 0
graz.data$pred.popn.24 <- graz.data$D1PredPop*(exp(graz.data$pred.mu.72)-1)/(graz.data$pred.mu.72) # mean predator population
graz.data$prey.i.24 <- graz.data$prey.g.72*graz.data$prey.popn.24/graz.data$pred.popn.24 # ingestion = mean prey * grazing rate / mean predator


graz.data$prey.i.24.mean <- graz.data$prey.g.72.mean*graz.data$prey.popn.24/graz.data$pred.popn.24 # ingestion = mean prey * grazing rate / mean predator

# Finally we'll manually fill in grazing and ingestion rates of 0 for 0 prey treatments.
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.g.72 <- 0
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.i.24 <- 0

graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.g.72.mean <- 0
graz.data[graz.data$TargetOchConc==0&graz.data$Pred!='Control',]$prey.i.24.mean <- 0

Ingestion vs population size figures

Let’s look at our results! A scatterplot lets us see our results as a function of the mean prey concentration. We can even work on fitting our data with a Holling Type II (saturating) functional response. Here, I’m also choosing to fit a line to the most linear part of the data to estimate the predator’s attack rate.

par(mar=c(2,2.5,1,1),mfrow=c(2,4))

linear.fit.maxima <- c(70000, 50000, 40000, 25000, 155000,90000, 60000, 80000)
saturating.fit.maxima <- c(100000, 100000, 90000, 135000, 155000,250000, 170000, 110000)


attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
  
  # Add labeling
  text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
  if(i == 1){
    legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
  }
  
  # Add xenic data
  subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
  
  # Add linear model for all data
  subdat1.trim <- subdat1[subdat1$prey.popn.24 < linear.fit.maxima[i],]
  lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = linear.fit.maxima[i])
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=1.5,col=col.strains[i])
  #abline(lm1)
  
  attack.ests$a.linear[i] <- summary(lm1)$parameters[1]
  attack.ests$a.linear.se[i] <- summary(lm1)$parameters[2]



  # Add saturating model for all data
  subdat1.trim <- subdat1[subdat1$prey.popn.24 < saturating.fit.maxima[i],]
  a.est <- summary(lm(subdat1.trim[1:7,]$prey.i.24~subdat1.trim[1:7,]$prey.popn.24))$coefficients[2,1] # estimate initial slope
  h.est <- 1/50 # estimate handling time

  fit1 <- nls(prey.i.24 ~ a*prey.popn.24/(1+a*h*prey.popn.24),data=subdat1.trim,start = list(a = a.est, h = h.est))
  a.fit <- summary(fit1)$parameters[1,1]
  h.fit <- summary(fit1)$parameters[2,1]

  xcoords <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out = 100)
  ycoords <- a.fit*xcoords/(1+a.fit*h.fit*xcoords)
  lines(xcoords,ycoords,lwd=2,col='black',lty=2)

  attack.ests$a.sat[i] <- summary(fit1)$parameters[1,1]
  attack.ests$a.sat.se[i] <- summary(fit1)$parameters[1,2]
  attack.ests$h.sat[i] <- summary(fit1)$parameters[2,1]
  attack.ests$h.sat.se[i] <- summary(fit1)$parameters[2,2]
}

Fitting linear models with cutoffs

par(mar=c(2,2.5,1,1),mfrow=c(2,4))

linear.fit.maxima <- c(70000, 50000, 40000, 25000, 155000,90000, 60000, 80000)
linear.fit.maxima.x <- c(150000, 100000, 50000, 25000, 155000,120000, 600000, 120000)
linear.fit.maxima.r <- c(70000, 100000, 40000, 150000, 155000,250000, 200000, 120000)
saturating.fit.maxima <- c(100000, 100000, 90000, 135000, 155000,250000, 170000, 110000)


attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.x.linear <- NaN
attack.ests$a.x.linear.se <- NaN
attack.ests$a.r.linear <- NaN
attack.ests$a.r.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
  
  # Add labeling
  text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
  if(i == 1){
    legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
  }
  
  # Add xenic data
  subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
  
  # Add linear model for all data
  subdat1.trim <- subdat1[subdat1$prey.popn.24 < linear.fit.maxima[i],]
  lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = linear.fit.maxima[i])
  yset <- summary(lm1)$parameters[1] * xset
  #lines(xset,yset,lwd=1.5,col=col.strains[i])
  #abline(lm1)
  
  attack.ests$a.linear[i] <- summary(lm1)$parameters[1]
  attack.ests$a.linear.se[i] <- summary(lm1)$parameters[2]
  
  
  # Add linear model for xenic data
  subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i] ,]
  subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i] & subdat1$prey.popn.24 < linear.fit.maxima.x[i],]
  lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=1.5,col=col.strains[i],lty=2)

  attack.ests$a.x.linear[i] <- summary(lm1)$parameters[1]
  attack.ests$a.x.linear.se[i] <- summary(lm1)$parameters[2]
  
  # Add linear model for rice data
  subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
  subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i] & subdat1$prey.popn.24 < linear.fit.maxima.r[i],]
  lm1 <- nls(prey.i.24 ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T))
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=1.5,col=col.strains[i],lty=1)

  attack.ests$a.r.linear[i] <- summary(lm1)$parameters[1]
  attack.ests$a.r.linear.se[i] <- summary(lm1)$parameters[2]



  # Add saturating model for all data
  subdat1.trim <- subdat1[subdat1$prey.popn.24 < saturating.fit.maxima[i],]
  a.est <- summary(lm(subdat1.trim[1:7,]$prey.i.24~subdat1.trim[1:7,]$prey.popn.24))$coefficients[2,1] # estimate initial slope
  h.est <- 1/50 # estimate handling time

  fit1 <- nls(prey.i.24 ~ a*prey.popn.24/(1+a*h*prey.popn.24),data=subdat1.trim,start = list(a = a.est, h = h.est))
  a.fit <- summary(fit1)$parameters[1,1]
  h.fit <- summary(fit1)$parameters[2,1]

  xcoords <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out = 100)
  ycoords <- a.fit*xcoords/(1+a.fit*h.fit*xcoords)
  #lines(xcoords,ycoords,lwd=2,col='black',lty=2)

  attack.ests$a.sat[i] <- summary(fit1)$parameters[1,1]
  attack.ests$a.sat.se[i] <- summary(fit1)$parameters[1,2]
  attack.ests$h.sat[i] <- summary(fit1)$parameters[2,1]
  attack.ests$h.sat.se[i] <- summary(fit1)$parameters[2,2]
}

Ingestion vs abundance, statistically robust

Adjust with breakpoint analysis to make cutoffs more statistically robust

par(mar=c(2,2.5,.7,.5))

rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
  for(j in 1:4){
    lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
  }
}
layout(lay.mat)

# Create a data frame to hold the results
attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.x.linear <- NaN
attack.ests$a.x.linear.se <- NaN
attack.ests$a.r.linear <- NaN
attack.ests$a.r.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
  
  # Add labeling
  text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
  if(i == 8){
    legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
  }
  if(i == 1){
    mtext(expression(paste('Ingestion rate (',italic('Ochromonas'),' · ',italic('O. marina')^-1,' · ',d^-1,')')),side=2,cex=.8,at=-6,line=2.5)
  }
  if(i == 6){
    mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
  }
  
  # Add xenic data
  subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
  
  # Add linear model for all data
  ## Test 3 types of model fits
  lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1) # +0 forces y-intercept through 0
  o1 <- segmented(lm1)
  o2 <- segmented(lm1,npsi=2)
  
  ## Select amongst them based on smallest p-value
  seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
  
  if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
    # Plot the best-fit model; skipping here
    # Record the data
    attack.ests$a.linear[i] <- summary(lm1)$coefficients[1,1]
    attack.ests$a.linear.se[i] <- summary(lm1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
    # Plot the best-fit model; skipping here
    # Record the data
    attack.ests$a.linear[i] <- summary(o1)$coefficients[1,1]
    attack.ests$a.linear.se[i] <- summary(o1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
    # Plot the best-fit model; skipping here
    # Record the data
    attack.ests$a.linear[i] <- summary(o2)$coefficients[1,1]
    attack.ests$a.linear.se[i] <- summary(o2)$coefficients[1,2]
  }

  # Add linear model for xenic data
  ## Test 3 types of model fits
  lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==xenic.strains[i] ,])
  o1 <- segmented(lm1)
  o2 <- segmented(lm1,npsi=2)
  
  ## Select amongst them based on smallest p-value
  seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
  
  if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
    # Plot the best-fit model
    xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==xenic.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
    ycoords <- xcoords*summary(lm1)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2,lty=2,col=col.strains[i])
    # Record the data
    attack.ests$a.x.linear[i] <- summary(lm1)$coefficients[1,1]
    attack.ests$a.x.linear.se[i] <- summary(lm1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
    # Plot the best-fit model
    plot.segmented(o1,add=TRUE,col=col.strains[i],lty=2,lwd=2)
    # Record the data
    attack.ests$a.x.linear[i] <- summary(o1)$coefficients[1,1]
    attack.ests$a.x.linear.se[i] <- summary(o1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
    # Plot the best-fit model
    plot.segmented(o2,add=TRUE,col=col.strains[i],lty=2,lwd=2)
    # Record the data
    attack.ests$a.x.linear[i] <- summary(o2)$coefficients[1,1]
    attack.ests$a.x.linear.se[i] <- summary(o2)$coefficients[1,2]
  }
  
  
  # Add linear model for rice data
  ## Test 3 types of model fits
  lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==rice.strains[i] ,])
  o1 <- segmented(lm1)
  if(i == 6){
    o2 <- o1
  } else{
  o2 <- segmented(lm1,npsi=2)}
  
  ## Select amongst them based on smallest p-value
  seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
  
  if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
    # Plot the best-fit model
    xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==rice.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
    ycoords <- xcoords*summary(lm1)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
    # Record the data
    attack.ests$a.r.linear[i] <- summary(lm1)$coefficients[1,1]
    attack.ests$a.r.linear.se[i] <- summary(lm1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
    # Plot the best-fit model
    plot.segmented(o1,add=TRUE,col=col.strains[i],lty=1,lwd=2)
    # Record the data
    attack.ests$a.r.linear[i] <- summary(o1)$coefficients[1,1]
    attack.ests$a.r.linear.se[i] <- summary(o1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
    # Plot the best-fit model
    plot.segmented(o2,add=TRUE,col=col.strains[i],lty=1,lwd=2)
    # Record the data
    attack.ests$a.r.linear[i] <- summary(o2)$coefficients[1,1]
    attack.ests$a.r.linear.se[i] <- summary(o2)$coefficients[1,2]
  }
}

Including just the first line segment

par(mar=c(2,2.5,.7,.5))

rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
  for(j in 1:4){
    lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
  }
}
layout(lay.mat)

# Create a data frame to hold the results
attack.ests <- as.data.frame(xenic.strains)
colnames(attack.ests) <-'Strain'
attack.ests$a.linear <- NaN
attack.ests$a.linear.se <- NaN
attack.ests$a.x.linear <- NaN
attack.ests$a.x.linear.se <- NaN
attack.ests$a.r.linear <- NaN
attack.ests$a.r.linear.se <- NaN
attack.ests$a.sat <- NaN
attack.ests$a.sat.se <- NaN
attack.ests$h.sat <- NaN
attack.ests$h.sat.se <- NaN


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(0,max(subdat1$prey.i.24,na.rm=T)),las=1)
  
  # Add labeling
  text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$prey.i.24,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
  if(i == 1){
    legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
    mtext(expression(paste('Ingestion rate (',italic('Ochromonas'),' · ',italic('O. marina')^-1,' · ',d^-1,')')),side=2,cex=.8,at=-6,line=2.5)
  }
  if(i == 6){
    mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
  }
  
  # Add xenic data
  subdat.x <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.x$prey.popn.24,subdat.x$prey.i.24,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,]
  points(subdat.r$prey.popn.24,subdat.r$prey.i.24, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
  
  # Add linear model for all data
  ## Test 3 types of model fits
  lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1) # +0 forces y-intercept through 0
  o1 <- segmented(lm1)
  o2 <- segmented(lm1,npsi=2)
  
  ## Select amongst them based on smallest p-value
  seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
  
  if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
    # Plot the best-fit model; skipping here
    # Record the data
    attack.ests$a.linear[i] <- summary(lm1)$coefficients[1,1]
    attack.ests$a.linear.se[i] <- summary(lm1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
    # Plot the best-fit model; skipping here
    # Record the data
    attack.ests$a.linear[i] <- summary(o1)$coefficients[1,1]
    attack.ests$a.linear.se[i] <- summary(o1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
    # Plot the best-fit model; skipping here
    # Record the data
    attack.ests$a.linear[i] <- summary(o2)$coefficients[1,1]
    attack.ests$a.linear.se[i] <- summary(o2)$coefficients[1,2]
  }

  # Add linear model for xenic data
  ## Test 3 types of model fits
  lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==xenic.strains[i] ,])
  o1 <- segmented(lm1)
  o2 <- segmented(lm1,npsi=2)
  
  ## Select amongst them based on smallest p-value
  seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
  
  if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
    # Plot the best-fit model
    xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==xenic.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
    ycoords <- xcoords*summary(lm1)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2.5,lty=2,col=col.strains[i])
    # Record the data
    attack.ests$a.x.linear[i] <- summary(lm1)$coefficients[1,1]
    attack.ests$a.x.linear.se[i] <- summary(lm1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
    # Plot the best-fit model
    #plot.segmented(o1,add=TRUE,col=col.strains[i],lty=2,lwd=2)
    xcoords <- seq(from = 0, to = o1$psi[1,2],length.out=100)
    ycoords <- xcoords*summary(o1)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2.5,lty=2,col=col.strains[i])
    # Record the data
    attack.ests$a.x.linear[i] <- summary(o1)$coefficients[1,1]
    attack.ests$a.x.linear.se[i] <- summary(o1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
    # Plot the best-fit model
    xcoords <- seq(from = 0, to = o2$psi[1,2],length.out=100)
    ycoords <- xcoords*summary(o2)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2.5,lty=2,col=col.strains[i])
    # Record the data
    attack.ests$a.x.linear[i] <- summary(o2)$coefficients[1,1]
    attack.ests$a.x.linear.se[i] <- summary(o2)$coefficients[1,2]
  }
  
  
  # Add linear model for rice data
  ## Test 3 types of model fits
  lm1 <- lm(prey.i.24 ~ prey.popn.24+0, data=subdat1[subdat1$Strain==rice.strains[i] ,])
  o1 <- segmented(lm1)
  if(i == 6){
    o2 <- o1
  } else{
  o2 <- segmented(lm1,npsi=2)}
  
  ## Select amongst them based on smallest p-value
  seg.pvals <- c(summary(lm1)$coefficients[1,4],summary(o1)$coefficients[1,4],summary(o2)$coefficients[1,4])
  
  if(which(seg.pvals==min(seg.pvals))[1]==1){ # If the standard linear model is the best
    # Plot the best-fit model
    xcoords <- seq(from = 0, to = max(subdat1[subdat1$Strain==rice.strains[i],]$prey.popn.24,na.rm=T),length.out=100)
    ycoords <- xcoords*summary(lm1)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
    # Record the data
    attack.ests$a.r.linear[i] <- summary(lm1)$coefficients[1,1]
    attack.ests$a.r.linear.se[i] <- summary(lm1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==2){ # If the 1 breakpoint model is the best
    # Plot the best-fit model
    #plot.segmented(o1,add=TRUE,col=col.strains[i],lty=2,lwd=2)
    xcoords <- seq(from = 0, to = o1$psi[1,2],length.out=100)
    ycoords <- xcoords*summary(o1)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
    # Record the data
    attack.ests$a.r.linear[i] <- summary(o1)$coefficients[1,1]
    attack.ests$a.r.linear.se[i] <- summary(o1)$coefficients[1,2]
  }
  if(which(seg.pvals==min(seg.pvals))[1]==3){ # If the 2 breakpoint model is the best
    # Plot the best-fit model
    xcoords <- seq(from = 0, to = o2$psi[1,2],length.out=100)
    ycoords <- xcoords*summary(o2)$coefficients[1,1]
    lines(xcoords,ycoords,lwd=2,lty=1,col=col.strains[i])
    # Record the data
    attack.ests$a.r.linear[i] <- summary(o2)$coefficients[1,1]
    attack.ests$a.r.linear.se[i] <- summary(o2)$coefficients[1,2]
  }
  
}

Barplot of attack rates of O. marina on Ochromonas

par(mar=c(4,4.2,1,1))
barplot(attack.ests$a.linear*1000,las=1, ylim=c(0,max(attack.ests$a.linear+attack.ests$a.linear.se))*1001,names=xenic.strains,xlab='Strain',ylab=expression(paste(italic('O. marina'),' clearance rate (μL ',d^-1,")")), border="black", col = col.strains.bg)
xcoords <- seq(from = 0.75, to = 9.1, length.out=8)
arrows(xcoords,attack.ests$a.linear*1000-attack.ests$a.linear.se*1000,xcoords, attack.ests$a.linear*1000+attack.ests$a.linear.se*1000, code = 3, angle = 90, length = 0.05, col= 'black')

subdat <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.r.linear,attack.ests$a.r.linear.se))
colnames(subdat) <- c('Strain','attack','se')
subdat$Bact <- 'R'
subdat2 <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.x.linear,attack.ests$a.x.linear.se))
colnames(subdat2) <- c('Strain','attack','se')
subdat2$Bact <- 'X'

subdat <- rbind(subdat2,subdat)
subdat <- subdat[order(subdat$Strain),]
barplot(subdat$attack*1000,las=1, ylim=c(0,max(subdat$attack+subdat$se))*1001,xlab='Strain',ylab=expression(paste(italic('O. marina'),' clearance rate (μL · ',d^-1,")")), border=rep(col.strains,each=2),lwd=1, col = mega.col.strains.white,space=c(0.7,.1))
xcoords1 <- seq(from = 1.2, to = 20.8, length.out=8)
xcoords2 <- seq(from = 2.3, to = 21.9, length.out = 8)
xcoords <- sort(c(xcoords1,xcoords2))
arrows(xcoords,subdat$attack*1000-subdat$se*1000,xcoords, subdat$attack*1000+subdat$se*1000, code = 3, angle = 90, length = 0.05, col= 'black')
axis(side=1, at = colMeans(rbind(xcoords1,xcoords2)), labels=xenic.strains, padj=-.5 )
legend(x='topright',inset=c(0,0),legend=c('Low Bacteria','High Bacteria'),pch=22,pt.cex=2,bty='n',pt.bg=c('white', 'gray10'), col=c('gray10','gray10'))

# Adding t-test results
barplot(subdat$attack*1000,las=1, ylim=c(0,max(subdat$attack+subdat$se))*1001,xlab='Strain',ylab=expression(paste(italic('O. marina'),' clearance rate (μL · ',d^-1,")")), border=rep(col.strains,each=2),lwd=1, col = mega.col.strains.white,space=c(0.7,.1))
xcoords1 <- seq(from = 1.2, to = 20.8, length.out=8)
xcoords2 <- seq(from = 2.3, to = 21.9, length.out = 8)
xcoords <- sort(c(xcoords1,xcoords2))
xcoords.mean <- colMeans(rbind(xcoords1,xcoords2))
arrows(xcoords,subdat$attack*1000-subdat$se*1000,xcoords, subdat$attack*1000+subdat$se*1000, code = 3, angle = 90, length = 0.05, col= 'black')
axis(side=1, at = xcoords.mean, labels=xenic.strains, padj=-.5 )
legend(x='topright',inset=c(0,0),legend=c('Low Bacteria','High Bacteria'),pch=22,pt.cex=2,bty='n',pt.bg=c('white', 'gray10'), col=c('gray10','gray10'))

pval.text <- rep('',8)
for(i in 1:length(xenic.strains)){
  subdat2 <- subdat[subdat$Strain==xenic.strains[i],]
  max.se <- max(subdat2$se)
  delta.attack <- max(subdat2$attack) - min(subdat2$attack)
  num.ses <- delta.attack / max.se
  if(num.ses > 3){
    pval.text[i] <- '*'
  }
  if(num.ses > 4){
    pval.text[i] <- '**'
  }
  if(num.ses > 5){
    pval.text[i] <- '***'
  }
}
ycoords <- tapply(subdat$attack*1000,subdat$Strain,FUN='max') + .04
text(xcoords.mean,ycoords,pval.text,cex=.75)

O. marina performance

relative growth rates

Raw growth rates

par(mar=c(2,2.5,.7,.5))

rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
  for(j in 1:4){
    lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
  }
}
layout(lay.mat)


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(min(subdat1$pred.mu.72,na.rm=T),max(subdat1$pred.mu.72,na.rm=T)),las=1)
  abline(h=0,lty=2)
  
  # Add xenic data
  subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
  points(subdat.x$prey.popn.24,subdat.x$pred.mu.72,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
  points(subdat.r$prey.popn.24,subdat.r$pred.mu.72, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
  
  # Add labeling
  text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$pred.mu.72,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
  if(i == 1){
    legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
  }
  if(i == 1){
    mtext(expression(paste(italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.6,line=2.5)
  }
  if(i == 6){
    mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
  }
}

Relativized growth rates

# First, calculate the mean growth rates of O. marina for each experiment
# Create a unique identifier for each experiment
graz.data$Strain.D1Date <- paste(graz.data$Strain,graz.data$D1Date,sep='.')

# Compute O. marina control (0 prey) growth rates for each experiment
ref.Omarina.growthrates <- summarySE(data=graz.data[graz.data$TargetOchConc==0 & !is.na(graz.data$pred.mu.72),],'pred.mu.72',groupvars='Strain.D1Date')

# Second, map the O. marina control growth rates back to their relevant experiment
graz.data$pred.mu.72.ctrl <- ref.Omarina.growthrates$pred.mu.72[match(graz.data$Strain.D1Date,ref.Omarina.growthrates$Strain.D1Date)]

# Third, compute the relativized growth rates
graz.data$pred.mu.72.rel <- graz.data$pred.mu.72-graz.data$pred.mu.72.ctrl

par(mar=c(2,2.5,.7,.5))

rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
  for(j in 1:4){
    lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
  }
}
layout(lay.mat)

 growth.ests <- as.data.frame(xenic.strains)
 colnames(growth.ests) <-'Strain'
 growth.ests$a.linear <- NaN
 growth.ests$a.linear.se <- NaN
 growth.ests$a.all.linear <- NaN
 growth.ests$a.all.linear.se <- NaN
 growth.ests$a.x.linear <- NaN
 growth.ests$a.x.linear.se <- NaN
 growth.ests$a.r.linear <- NaN
 growth.ests$a.r.linear.se <- NaN
 growth.ests$a.sat <- NaN
 growth.ests$a.sat.se <- NaN
 growth.ests$h.sat <- NaN
 growth.ests$h.sat.se <- NaN


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])
  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.popn.24,na.rm=T)), ylim=c(min(subdat1$pred.mu.72.rel,na.rm=T),max(subdat1$pred.mu.72.rel,na.rm=T)),las=1)
  abline(h=0,lty=1,col='gray50')
  
  # Add labeling
  if(i %in% c(1,2,5,6,7,8)){
    text(x = 0-.06*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)
    }
  if(i %in% c(3,4)){
    text(x = 1.05*max(subdat1$prey.popn.24,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 2)    
  }
  if(i == 8){
    legend(x='topright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2)
  }
  if(i == 1){
    mtext(expression(paste('Relativized ', italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.2,line=2.5)
  }
  if(i == 6){
    mtext(expression(paste('Prey availability (',italic('Ochromonas'),' · ',mL^-1, ')')),side=1,line=2.7,cex=.8,at=275000)
  }
  
  # Add xenic data
  subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
  points(subdat.x$prey.popn.24,subdat.x$pred.mu.72.rel,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
  points(subdat.r$prey.popn.24,subdat.r$pred.mu.72.rel, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])
  
  # Add linear model for subset of data
  subdat1.trim <- subdat1[subdat1$prey.popn.24 < linear.fit.maxima[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = linear.fit.maxima[i])
  yset <- summary(lm1)$parameters[1] * xset
  #lines(xset,yset,lwd=1.5,col=col.strains[i])

  growth.ests$a.linear[i] <- summary(lm1)$parameters[1]
  growth.ests$a.linear.se[i] <- summary(lm1)$parameters[2]
  
  
  # Add linear model for all data
  subdat1.trim <- subdat1
  lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  #lines(xset,yset,lwd=1.5,col='black',lty=2)

  growth.ests$a.all.linear[i] <- summary(lm1)$parameters[1]
  growth.ests$a.all.linear.se[i] <- summary(lm1)$parameters[2]


  # Add linear model for xenic data
  subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=2.5,col=col.strains[i],lty=2)

  growth.ests$a.x.linear[i] <- summary(lm1)$parameters[1]
  growth.ests$a.x.linear.se[i] <- summary(lm1)$parameters[2]
  
  # Add linear model for rice data
  subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.popn.24, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.popn.24,na.rm=T))
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=2,col=col.strains[i],lty=1)

  growth.ests$a.r.linear[i] <- summary(lm1)$parameters[1]
  growth.ests$a.r.linear.se[i] <- summary(lm1)$parameters[2]

}

growth vs ingestion

graz.data$Strain.numeric <- as.numeric(str_extract_all(graz.data$Strain, "\\d+"))
graz.data$Bact <- 'X'
graz.data[graz.data$Strain%in%c('584R','590R','1148R','1150R','1391R','1392R','1393R','2951R'),]$Bact <- "R"

graz.data$Strain.Bact <- paste(graz.data$Strain.numeric,graz.data$Bact,sep='.')
Barbaglia.data$Strain.Bact <- paste(Barbaglia.data$Strain,Barbaglia.data$Bact.short,sep='.')

# Append information for Ochromonas carbon content in pg C per cell
graz.data$C.perOch <- Barbaglia.data$C.perOch[match(graz.data$Strain.Bact,Barbaglia.data$Strain.Bact)]

# Convert ingestion rate from cells to C, divide by 1000 to get ng per day
graz.data$prey.i.24.C <- graz.data$prey.i.24*graz.data$C.perOch/1000

# Append information for Ochromonas carbon content in pg N per cell
graz.data$N.perOch <- Barbaglia.data$N.perOch[match(graz.data$Strain.Bact,Barbaglia.data$Strain.Bact)]

# Convert ingestion rate from cells to C, divide by 1000 to get ng per day
graz.data$prey.i.24.N <- graz.data$prey.i.24*graz.data$N.perOch/1000

par(mar=c(2,2.5,.7,.5))

rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
  for(j in 1:4){
    lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
  }
}
layout(lay.mat)


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])

  
  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.i.24.C,na.rm=T)), ylim=c(min(subdat1$pred.mu.72.rel,na.rm=T),max(subdat1$pred.mu.72.rel,na.rm=T)),las=1)
  abline(h=0,lty=1,col='gray50')
  
  # Add xenic data
  subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
  points(subdat.x$prey.i.24.C,subdat.x$pred.mu.72.rel,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
  points(subdat.r$prey.i.24.C,subdat.r$pred.mu.72.rel, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])

  # Add linear model for xenic data
  subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.C, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.C,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=2.5,col=col.strains[i],lty=2)

  # Add linear model for rice data
  subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.C, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.C,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=2,col=col.strains[i],lty=1)
  
  # Add labeling
  if(i %in% c(1,5:8)){
    text(x = 0-.06*max(subdat1$prey.i.24.C,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)}
  if(i %in% c(2:4)){
    text(x = 1.06*max(subdat1$prey.i.24.C,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 2)    
  }
  if(i == 8){
    legend(x='bottomright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2, lty = c(2,1))
  }
  if(i == 1){
    mtext(expression(paste('Relativized ', italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.18,line=2.5)
  }
  if(i == 6){
    mtext(expression(paste('Ingestion rate (ng ',italic('Ochromonas'),' C · ',italic('O. marina')^-1,' · ',d^-1,')')),side=1,line=2.7,cex=.8,at=4.5)
  }

}

par(mar=c(2,2.5,.7,.5))

rowsperpanel <- 5
colsperpanel <- 5
lay.mat <- matrix(data=NA,nrow=2*rowsperpanel+1,ncol=4*colsperpanel+1)
lay.mat[,1] <- 9
lay.mat[2*rowsperpanel+1,] <- 10
for(i in 1:2){
  for(j in 1:4){
    lay.mat[(1+(i-1)*rowsperpanel):(i*rowsperpanel),(2+(j-1)*colsperpanel):(1+j*colsperpanel)] <- 4*(i-1) + j
  }
}
layout(lay.mat)


for(i in 1:length(xenic.strains)){
  subdat1 <- rbind(graz.data[graz.data$TargetOchConc!='Control'& graz.data$Strain==xenic.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,],graz.data[graz.data$Pred!='Control'& graz.data$Strain==rice.strains[i] & graz.data$OldDat2!='Y' & graz.data$prey.i.24 >= 0,])

  # Make the plot
  plot(1,1,type='n',xlab='',ylab='',xlim=c(0,max(subdat1$prey.i.24.N,na.rm=T)), ylim=c(min(subdat1$pred.mu.72.rel,na.rm=T),max(subdat1$pred.mu.72.rel,na.rm=T)),las=1)
  abline(h=0,lty=1,col='black')
  
  # Add xenic data
  subdat.x <- subdat1[subdat1$Strain==xenic.strains[i] ,]
  points(subdat.x$prey.i.24.N,subdat.x$pred.mu.72.rel,pch=21, cex=1.6, col=col.strains[i], bg='white')# bg= col.strains.bg[i])
  
  # Add rice data
  subdat.r <- subdat1[subdat1$Strain==rice.strains[i],]
  points(subdat.r$prey.i.24.N,subdat.r$pred.mu.72.rel, pch=24, cex=1.2, col = col.strains[i], bg = col.strains[i])

  # Add linear model for xenic data
  subdat1.trim <- subdat1[subdat1$Strain==xenic.strains[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.N, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.N,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=2.5,col=col.strains[i],lty=2)

  # Add linear model for rice data
  subdat1.trim <- subdat1[subdat1$Strain==rice.strains[i],]
  lm1 <- nls(pred.mu.72.rel ~ a*prey.i.24.N, data=subdat1.trim, start = list(a = 0.001))
  xset <- seq(from = 0, to = max(subdat1.trim$prey.i.24.N,na.rm=T),length.out=100)
  yset <- summary(lm1)$parameters[1] * xset
  lines(xset,yset,lwd=2,col=col.strains[i],lty=1)


  # Add labeling
  if(i %in% c(1,5,6,8)){
    text(x = 0-.06*max(subdat1$prey.i.24.N,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 4)}
  if(i %in% c(2:4,7)){
    text(x = 1.06*max(subdat1$prey.i.24.N,na.rm=T), y = max(subdat1$pred.mu.72.rel,na.rm=T)*0.97, xenic.strains[i], col = col.strains[i],cex=1.3, pos = 2)    
  }
  if(i == 8){
    legend(x='bottomright', inset=c(0.02,0.03), legend=c('Low','High'), pch=c(21,24), col=col.strains[i], pt.bg=c('white',col.strains[i]), bty='y', pt.cex=1.6,1.2, lty = c(2,1))
  }
  if(i == 1){
    mtext(expression(paste('Relativized ', italic('O. marina'),' growth rate (',d^-1,')')),side=2,cex=.8,at=-.18,line=2.5)
  }
  if(i == 6){
    mtext(expression(paste('Ingestion rate (ng ',italic('Ochromonas'),' N · ',italic('O. marina')^-1,' · ',d^-1,')')),side=1,line=2.7,cex=.8,at=.9)
  }

  
}

THE BIG RESULT: testing for correlation between Ochromonas palatability and Barbaglia et al. traits

Barbaglia.data$Bact.short <- factor(Barbaglia.data$Bact.short,levels=c('X','R'))


# Useful features:
## C.from.grazing.perCperday = carbon obtained from grazing per Ochromonas carbon per day
## P_I.perC = carbon obtained from photosynthesis per Ochromonas carbon per day

Grazing

#Barbaglia.data$Strain.Bact.Light
par(mar=c(4,4,1,1))

attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'Best.a', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab='Ochromonas attack rate', ylab='O. marina attack rate',type='n',ylim=c(0.3,1.6),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=c(24,21)[as.factor(attack.summ$Bact.short)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' clearance rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')

Chlorophyll

par(mar=c(4,4,1,1))

attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'chl.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' chlorophyll per cell')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.8),xlim=c(0.12,0.34))
arrows(attack.summ$chl.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$chl.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$chl.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')

Both together

par(mar=c(4,4,1,0.5),mfrow=c(1,2))

attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'Best.a', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' clearance rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')



par(mar=c(4,2.5,1,2))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'chl.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' chlorophyll per cell')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95)) #,xlim=c(0.12,0.34)
arrows(attack.summ$chl.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$chl.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$chl.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$chl.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)

#legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')

Total carbon yield

par(mar=c(4,4,1,1))

## C.from.grazing.perCperday = carbon obtained from grazing per Ochromonas carbon per day
## P_I.perC = carbon obtained from photosynthesis per Ochromonas carbon per day

Barbaglia.data$TotC.perC <- Barbaglia.data$C.from.grazing.perCperday + Barbaglia.data$P_I.perC
Barbaglia.data$TotC.percell <- Barbaglia.data$TotC.perC * Barbaglia.data$C.perOch

attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'TotC.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' carbon yield rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0,250))
arrows(attack.summ$TotC.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$TotC.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$TotC.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')

Combining this with the attack rate figure

par(mar=c(4,4,1,0.5),mfrow=c(1,2))

attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'Best.a', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' clearance rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0.3e-6,4.7e-6))
arrows(attack.summ$Best.a-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$Best.a+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$Best.a,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$Best.a,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$Best.a,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)
legend(x='topright',inset=c(0.01,0.01),legend=xenic.strains,pch=pch.strains, cex=.8, pt.cex=cex.strains, col=col.strains, pt.bg=col.strains, title='Strain')



par(mar=c(4,2.5,1,2))
attack.summ <- summarySE(data=Barbaglia.data[Barbaglia.data$Light==100,],'TotC.percell', groupvars=c('Strain','Bact.short'))
attack.summ$O.marina.a <- attack.ests$a.linear[match(attack.summ$Strain,attack.ests$Strain)]
attack.summ$O.marina.a.bact<-NaN
attack.summ$O.marina.a.bact.se<-NaN
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact <- attack.ests$a.x.linear[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='X',]$O.marina.a.bact.se <- attack.ests$a.x.linear.se[match(attack.summ[attack.summ$Bact.short=='X',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact <- attack.ests$a.r.linear[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]
attack.summ[attack.summ$Bact.short=='R',]$O.marina.a.bact.se <- attack.ests$a.r.linear.se[match(attack.summ[attack.summ$Bact.short=='R',]$Strain,attack.ests$Strain)]

pch.strains <- c(21,22,23,24,25,21,22,23)
cex.strains <- c(1.2,1.2,1.2,1,1,1.2,1.2,1.2)
plot(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,las=1, xlab=expression(paste(italic('Ochromonas'),' carbon yield rate')), ylab=expression(paste(italic('O. marina'),' clearance rate')),type='n',ylim=c(0.2,1.95),xlim=c(0,240))
arrows(attack.summ$TotC.percell-attack.summ$se,attack.summ$O.marina.a.bact*1000, attack.summ$TotC.percell+attack.summ$se,attack.summ$O.marina.a.bact*1000, code = 3, angle = 90, length = 0.03)
arrows(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 - attack.summ$O.marina.a.bact.se*1000, attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000 + attack.summ$O.marina.a.bact.se*1000, code = 3, angle = 90, length = 0.03)
for(i in 1:length(strains)){
  subdat <- attack.summ[attack.summ$Strain==strains[i],]
  lines(subdat$TotC.percell,subdat$O.marina.a.bact*1000,lty=2,col='gray50')
}
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains.white)
points(attack.summ$TotC.percell,attack.summ$O.marina.a.bact*1000,pch=pch.strains[as.factor(attack.summ$Strain)],cex = cex.strains[as.factor(attack.summ$Strain)], col = col.strains[as.factor(attack.summ$Strain)], bg= mega.col.strains)

T-tests

# Make a holding table
t.test.results <- as.data.frame(xenic.strains)
colnames(t.test.results) <- "Strain"
t.test.results$Och.graz.tval <- NaN
t.test.results$Och.graz.pval <- NaN
t.test.results$Oxy.graz.numSEs <- NaN

for(i in 1:length(xenic.strains)){ # Working strain by strain
strainchoice <- xenic.strains[i]

# Gina's data: Is there a difference in clearance rate between xenic and rice data?
# Negative t-value for Gina's data means the rice strains had *higher* clearance rates than the xenic ones
gina.subdat <- Barbaglia.data[Barbaglia.data$Light==100 & Barbaglia.data$Strain==strainchoice,]
t.test.results$Och.graz.tval[i] <- t.test(Best.a ~ Bact.short,data=gina.subdat)$statistic
t.test.results$Och.graz.pval[i] <- t.test(Best.a ~ Bact.short,data=gina.subdat)$p.value

}

# AG's data: Is there a difference in clearance rate between xenic and rice data?
# Negative test value means rice strains had *higher* clearance rates
subdat <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.r.linear,attack.ests$a.r.linear.se))
colnames(subdat) <- c('Strain','attack','se')
subdat$Bact <- 'R'
subdat2 <- as.data.frame(cbind(attack.ests$Strain,attack.ests$a.x.linear,attack.ests$a.x.linear.se))
colnames(subdat2) <- c('Strain','attack','se')
subdat2$Bact <- 'X'

subdat <- rbind(subdat2,subdat)
subdat <- subdat[order(subdat$Strain),]

for(i in 1:length(xenic.strains)){
  subdat2 <- subdat[subdat$Strain==xenic.strains[i],]
  max.se <- max(subdat2$se)
  delta.attack <- subdat2[subdat2$Bact=='X',]$attack - subdat2[subdat2$Bact=='R',]$attack
  t.test.results$Oxy.graz.numSEs[i] <- delta.attack / max.se
}

print(t.test.results)

##   Strain Och.graz.tval Och.graz.pval Oxy.graz.numSEs
## 1    584     3.0993720  0.0647563383       -8.778937
## 2    590    -0.2824363  0.7981552021       -2.951987
## 3   1148    -0.1434081  0.8937054002        1.835591
## 4   1150     0.6158596  0.5903486015        1.268617
## 5   1391    11.5989292  0.0003357685      -10.426857
## 6   1392     5.4398070  0.0075964682       -5.033688
## 7   1393     9.8461959  0.0006750262       -3.972854
## 8   2951    18.9398352  0.0009854816      -15.845488